Size Analysis of the Northeast US Finfish Community
Published
April 26, 2024
Modeling Change in Median Size and Community Spectra
The following tabs detail how I’m thinking about whether to use SST in Celsius or as anomalies to test our hypotheses about temperature’s impact on size and size structure:
SST Change Over Time
Landings Change Over Time
Landings Verification
Bart raised that his landings looked different, vdrify that we didn’t mess them up in data prep:
Landings and Anomalies Colinearity
Landings Clustering/Bi-modality
Spectra Slope Models
Contemporary SST/Landings Model
effect
group
term
estimate
std.error
statistic
fixed
NA
(Intercept)
-0.921360717
0.009839746
-93.6366360
fixed
NA
landings
0.018576979
0.010705595
1.7352589
fixed
NA
regionGeorges_Bank
-0.012148242
0.012621505
-0.9625034
fixed
NA
regionMid_Atlantic_Bight
-0.114884622
0.012637417
-9.0908308
fixed
NA
regionSouthern_New_England
-0.117012827
0.012378032
-9.4532658
fixed
NA
sst_anom
-0.004065130
0.009320448
-0.4361518
fixed
NA
landings:regionGeorges_Bank
-0.006458479
0.012945264
-0.4989067
fixed
NA
landings:regionMid_Atlantic_Bight
-0.031802414
0.013679938
-2.3247484
fixed
NA
landings:regionSouthern_New_England
-0.005036802
0.012402315
-0.4061179
ran_pars
year
sd__(Intercept)
0.022544228
NA
NA
ran_pars
Residual
sd__Observation
0.050777216
NA
NA
Significance
npar
AIC
LRT
Pr(Chi)
NA
-436.8383
NA
NA
sst_anom
1
-438.6385
0.1998047
0.6548785
landings:region
3
-435.7297
7.1086049
0.0685156
Fit & Predictions
Collinearity
SST Effect
Landings Effect
5-Year Smoothed Predictor Model
effect
group
term
estimate
std.error
statistic
fixed
NA
(Intercept)
-0.925286340
0.009426778
-98.15509603
fixed
NA
land_5
0.019084432
0.009964401
1.91526141
fixed
NA
regionGeorges_Bank
-0.008202934
0.011842074
-0.69269404
fixed
NA
regionMid_Atlantic_Bight
-0.111854877
0.011990018
-9.32899978
fixed
NA
regionSouthern_New_England
-0.116485618
0.011704724
-9.95201773
fixed
NA
anom_5
-0.001404379
0.014197197
-0.09891946
fixed
NA
land_5:regionGeorges_Bank
0.005290628
0.012412550
0.42623218
fixed
NA
land_5:regionMid_Atlantic_Bight
-0.033785136
0.013692523
-2.46741497
fixed
NA
land_5:regionSouthern_New_England
-0.011822140
0.012133137
-0.97436792
ran_pars
year
sd__(Intercept)
0.022792284
NA
NA
ran_pars
Residual
sd__Observation
0.050168060
NA
NA
Significance
npar
AIC
LRT
Pr(Chi)
NA
-442.7303
NA
NA
anom_5
1
-444.7201
0.0101711
0.9196679
land_5:region
3
-439.2276
9.5026805
0.0233029
Fit & Predictions
Collinearity
SST Effect
Landings Effect
Median Body Weight Models
How has median weight changed
Contemporary SST/Landings Model
effect
group
term
estimate
std.error
statistic
fixed
NA
(Intercept)
0.176728059
0.009359083
18.883053
fixed
NA
landings
0.021297703
0.010155283
2.097204
fixed
NA
regionGeorges_Bank
-0.019398093
0.011926868
-1.626420
fixed
NA
regionMid_Atlantic_Bight
-0.124714909
0.011948834
-10.437413
fixed
NA
regionSouthern_New_England
-0.121706242
0.011696896
-10.405003
fixed
NA
sst_anom
-0.002809302
0.008914145
-0.315151
fixed
NA
landings:regionGeorges_Bank
-0.013659094
0.012239620
-1.115974
fixed
NA
landings:regionMid_Atlantic_Bight
-0.020025091
0.012965115
-1.544536
fixed
NA
landings:regionSouthern_New_England
-0.018056535
0.011724685
-1.540044
ran_pars
year
sd__(Intercept)
0.022163153
NA
NA
ran_pars
Residual
sd__Observation
0.047978561
NA
NA
Fit & Predictions
Significance
npar
AIC
LRT
Pr(Chi)
NA
-452.6072
NA
NA
sst_anom
1
-454.5014
0.1058008
0.7449765
landings:region
3
-455.4081
3.1991481
0.3619278
Collinearity
SST Effect
Landings Effect
5-Year Smoothed Predictor Model
effect
group
term
estimate
std.error
statistic
fixed
NA
(Intercept)
0.17088659
0.008870858
19.2638172
fixed
NA
land_5
0.02683093
0.009353268
2.8686152
fixed
NA
regionGeorges_Bank
-0.01466959
0.011082561
-1.3236639
fixed
NA
regionMid_Atlantic_Bight
-0.12022181
0.011224651
-10.7105166
fixed
NA
regionSouthern_New_England
-0.11868733
0.010954052
-10.8350156
fixed
NA
anom_5
0.01107564
0.013407216
0.8260954
fixed
NA
land_5:regionGeorges_Bank
-0.00506618
0.011619881
-0.4359924
fixed
NA
land_5:regionMid_Atlantic_Bight
-0.02535527
0.012852623
-1.9727704
fixed
NA
land_5:regionSouthern_New_England
-0.02201549
0.011357699
-1.9383757
ran_pars
year
sd__(Intercept)
0.02197317
NA
NA
ran_pars
Residual
sd__Observation
0.04694883
NA
NA
Fit & Predictions
Significance
npar
AIC
LRT
Pr(Chi)
NA
-461.8534
NA
NA
anom_5
1
-463.1299
0.7234317
0.3950205
land_5:region
3
-461.3420
6.5113491
0.0892160
Collinearity
SST Effect
Landings Effect
Region Random Effects
What if we swap Random effects
I don’t really care about regions, so what if we used region as a random effect. It should also eat variability similar to year as a random effect and maybe help resolve the temperature or landings effects.
Just to see:
Contemporary SST/Landings Model
effect
group
term
estimate
std.error
df
statistic
p.value
fixed
NA
(Intercept)
0.109362822
0.028243777
145
3.8721033
0.0001627259
fixed
NA
landings
0.004371285
0.005859928
145
0.7459624
0.4568977897
fixed
NA
sst_anom
-0.003793204
0.008033209
145
-0.4721904
0.6375007698
ran_pars
survey_area
sd_(Intercept)
0.054752813
NA
NA
NA
NA
ran_pars
Residual
sd_Observation
0.053207811
NA
NA
NA
NA
Fit & Predictions
Significance
Df
AIC
LRT
Pr(>Chi)
NA
-440.8509
NA
NA
landings
1
-442.2872
0.5636696
0.4527855
sst_anom
1
-442.6271
0.2237871
0.6361694
Collinearity
SST Effect
Landings Effect
5-Year Smoothed Predictor Model
effect
group
term
estimate
std.error
statistic
fixed
NA
(Intercept)
-0.98269973
0.031874118
-30.8306483
fixed
NA
land_5
0.00173992
0.006155699
0.2826518
fixed
NA
anom_5
0.01589537
0.016403943
0.9689969
fixed
NA
land_5:anom_5
0.02932610
0.013546341
2.1648727
ran_pars
survey_area
sd__(Intercept)
0.06289653
NA
NA
ran_pars
Residual
sd__Observation
0.05544279
NA
NA
Significance
npar
AIC
LRT
Pr(Chi)
NA
-424.4753
NA
NA
land_5:anom_5
1
-421.7740
4.701356
0.0301388
Fit & Predictions
Collinearity
SST Effect
Landings Effect
Variance Partitioning
Source Code
---title: "Spectra and Median Size Mods"description: | Size Analysis of the Northeast US Finfish Communitydate: "`r Sys.Date()`"format: html: self-contained: true code-fold: true code-tools: true df-print: kable toc: true toc-depth: 2execute: echo: false warning: false message: false fig.align: center comment: ""---# Modeling Change in Median Size and Community Spectra```{r}library(gt)library(tidyverse)library(gmRi)library(tidyverse)library(ggeffects)library(zoo)library(lme4)library(patchwork)library(emmeans)library(rstatix)library(lmerTest)library(broom)library(broom.mixed)library(ggpmisc)library(performance)library(nlme)# Set themetheme_set(theme_gmri() +theme(plot.margin =margin(), legend.position ="bottom", legend.direction ="horizontal") )# levels for faceting areasarea_levels <-c("Northeast Shelf", "GoM", "GB", "SNE", "MAB")area_levels_long <-c("Northeast Shelf", "Gulf of Maine", "Georges Bank", "Southern New England", "Mid-Atlantic Bight")#--------------------## Set up#--------------------# Load the raw data df <-read_csv(here::here("data/size_and_spectra_model_data.csv"))# Get raw landingslandings_raw <-read_csv(here::here("data/unscaled_spectra_predictor_dataset_wide.csv")) %>%select(year, survey_area, landings_raw = landings)# Get the raw sst from oisstsst_raw =read_csv(here::here("data/unscaled_regional_sst.csv"))#------------# Put it all together and clarify the column namesdf <- df %>%select(-sst) %>%left_join(landings_raw) %>%left_join(sst_raw) %>%mutate(survey_area =factor( survey_area, levels = area_levels_long),region =str_replace_all(survey_area, "-| ", "_"),region =factor( region, levels =c("Gulf_of_Maine", "Georges_Bank", "Mid_Atlantic_Bight", "Southern_New_England")))# format year and anomaliesdf <- df %>%mutate(yr_num =as.numeric(as.character(year)),year =factor(year)) %>%group_by(survey_area) %>%arrange(survey_area, year) %>%mutate(anom_5 = zoo::rollapply( sst_anom, 5, mean, na.rm = T, align ="right", fill =NA),land_5_raw = zoo::rollapply( landings_raw, 5, mean, na.rm = T, align ="right", fill =NA))```The following tabs detail how I'm thinking about whether to use SST in Celsius or as anomalies to test our hypotheses about temperature's impact on size and size structure:### SST Change Over Time```{r}# Plot: How does temperature change over timesst_1 <- df %>%mutate(year =as.numeric(as.character(year))) %>%drop_na(sst) %>%ggplot(aes(year, sst, color = region)) +geom_point(size =1, alpha =0.6) +geom_smooth(method ="lm", linewidth =1, se = F) +scale_color_gmri() +labs(title ="Regional SST Increasing")# Plot: Anomaly Changessst_2 <- df %>%mutate(year =as.numeric(as.character(year))) %>%drop_na(sst_anom) %>%ggplot(aes(year, sst_anom, color = region)) +geom_point(size =1, alpha =0.6) +geom_smooth(method ="lm", linewidth =1, se = F) +scale_color_gmri() +labs(title ="Regional SST Anomalies Increasing")sst_1 / sst_2```### Landings Change Over Time```{r}# Plot: How does landings change over timedf %>%mutate(year =as.numeric(as.character(year))) %>%drop_na(landings_raw) %>%ggplot(aes(year, landings_raw, color = region)) +geom_point() +geom_smooth(method ="lm") +scale_color_gmri() +scale_y_continuous(labels = scales::label_comma()) +labs(title ="Regional Landings Declining",subtitle ="Mid-Atlantic Exception, primarily forage fish fishery (menhaden)")# Trying to keep landings: binned/categorical label```### Landings VerificationBart raised that his landings looked different, vdrify that we didn't mess them up in data prep:```{r}#| label: garfo-data-prep#| eval: true#### GARFO Landings ##### Load the GARFO landings data:# Landings of finfish* sheet 5res_path <-cs_path("res")landings <- readxl::read_xlsx(path =str_c(res_path, "GARFO_landings/KMills_landings by area 1964-2021_JUN 2022.xlsx"), sheet =5) %>%rename_all(tolower)# # Why no GB Landings in 2011?# landings %>% filter(# `stat area` %in% fish_zones$"Georges Bank", year == 2011)# Make a list of zones to roughly match the survey areas:fish_zones <-list("Gulf of Maine"=c(511:515, 464, 465),"Georges Bank"=c(521, 522, 525, 561, 562),"Southern New England"=c(611, 612, 613, 616, 526, 537, 538, 539),"Mid-Atlantic Bight"=c(614:615, 621, 622, 625, 626, 631, 632))# Add the labels into the landings data and remove what we don't need there:landings <- landings %>%mutate(survey_area =case_when(`stat area`%in% fish_zones$"Gulf of Maine"~"Gulf of Maine",`stat area`%in% fish_zones$"Georges Bank"~"Georges Bank",`stat area`%in% fish_zones$"Southern New England"~"Southern New England",`stat area`%in% fish_zones$"Mid-Atlantic Bight"~"Mid-Atlantic Bight")) %>%filter(survey_area %in%c("Georges Bank", "Gulf of Maine", "Southern New England", "Mid-Atlantic Bight")) %>%mutate(survey_area =factor(survey_area, area_levels_long))# Get Summarieslandings_summ <- landings %>%rename("weight_lb"=`landed lbs`,"live_lb"=`live lbs`) %>%group_by(year, survey_area) %>%summarise( across(.cols =c(value, weight_lb, live_lb), .fns =list(mean =~mean(.x , na.rm = T), total =~sum(.x , na.rm = T)), .names ="{.fn}_{.col}"), .groups ="drop") # Scale the landings to create an index by arealandings_summ <- landings_summ %>%group_by(survey_area) %>%mutate(total_wt_z =as.numeric(base::scale(total_weight_lb))) %>%ungroup()landings_summ %>%filter(year >=1970) %>%ggplot(aes(year, total_weight_lb, color = survey_area)) +geom_point() +geom_smooth(method ="lm") +scale_color_gmri() +scale_y_continuous(labels = scales::label_comma()) +labs(title ="Regional Landings Declining",subtitle ="Mid-Atlantic Exception, primarily forage fish fishery (menhaden)")```### Landings and Anomalies Colinearity```{r}# See roll anom and roll landings relationshipggplot(df, aes(anom_5, land_5_raw)) +geom_point() +facet_wrap(~region) +stat_poly_eq()# cor.test(df$sst_anom, df$landings_raw) #-.24# cor.test(df$anom_5, df$land_5_raw) # -.35# cor.test(df$sst_anom, df$landings) # -.46# cor.test(df$anom_5, df$land_5) # -.6```### Landings Clustering/Bi-modality```{r}# Exploring bimodality of landingsggplot(df, aes(landings)) +geom_density(aes(color = survey_area)) ``````{r}#| label: landings-clustering#| eval: false# clustering landings# for 2 clustersregion_mods <- df %>%drop_na(landings) %>%split(.$region) %>%imap(function(x,y){ mod =kmeans(drop_na(x,landings)$landings, centers =2)data.frame(region = y,year = x$year,cluster = mod$cluster )})# use betweens to assess which value for centers maximizes between group SS# evaluate n-clusters clust_range <-c(2:9)clust_num_mods <- clust_range %>%map(function(x){ landings_df <-drop_na(df,landings) mod =kmeans(landings_df$landings, centers = x)list("cluster_dat"=data.frame(year = landings_df$year,region = landings_df$region,cluster = mod$cluster),"betweens"=tibble("Between Group SS"= mod$betweens)) })# Pull performancemap_dfr(set_names(clust_num_mods, clust_range), ~pluck(.x, "betweens"), .id ="num_clusters") %>%ggplot(aes(num_clusters, `Between Group SS`)) +geom_point() +geom_line()# Pull the cluster info, make sure that each region has good representationdf %>%ggplot(aes(landings)) +geom_density() +facet_wrap(~region) +labs(title ="Landings (scaled) PDF")```---# Spectra Slope Models## Contemporary SST/Landings Model```{r}## 5-year smooth on anomalies## temperature impacts on growth and/or migration at large scales not instantaneous and integrate exposure over long time periods## # new model# From low to high correlation between landings and sst anomalies# Landings in raw form should be scaled# new_mod <- lme4::lmer(b ~ landings_raw * region + sst_anom + (1 | year), data = df) # new_mod <- lme4::lmer(b ~ land_5_raw * region + anom_5 + (1 | year), data = df)# To smooth or not to smooth, I vote nonew_mod <- lme4::lmer(b ~ landings * region + sst_anom + (1| year), data =drop_na(df, landings, sst_anom))# new_mod <- lme4::lmer(b ~ land_5 * region + anom_5 + (1 | year), data = df)# new_mod <- lme4::lmer(b ~ landings * region * sst_anom + (1 | year), data = df)# # Autocorrelative Error - Performs badly here# Ben bolker to thatnk, of course:# https://bbolker.github.io/mixedmodels-misc/notes/corr_braindump.html# new_mod <- nlme::lme(b ~ landings * region + sst_anom,# random = ~ 1 | year, # data = drop_na(df, landings, sst_anom),# correlation=corAR1(), method = "ML")# Summarytidy(new_mod) %>%gt()# # % variance of random effect# # year (intercept) / year (intercept) + residual# 0.0005082/(0.0005082 + 0.0025783) ```### Significance```{r}drop1(new_mod, test ="Chisq")```### Fit & Predictions```{r}qqnorm(resid(new_mod))qqline(resid(new_mod))ggplot(df, aes(x = yr_num, y = b, color = region)) +geom_point(alpha =0.5) +theme_classic() +scale_color_gmri() +geom_line(data =cbind(drop_na(df, landings, sst_anom), pred =predict(new_mod)), aes(y = pred), size =1) +# adding predicted line from mixed model theme(legend.position ="none",panel.spacing =unit(2, "lines")) # adding space between panels```### Collinearity```{r}# multicollinearity checkcol_check <-check_collinearity(new_mod)plot(col_check)```### SST Effect```{r}# SST Anomaliezsplot(ggpredict(new_mod, ~sst_anom+region), add.data = T, jitter =0.1)```### Landings Effect```{r}plot(ggpredict(new_mod, ~landings+region), add.data = T, jitter =0.1)```## 5-Year Smoothed Predictor Model```{r}# Both sst & landings - clears collinearity testlag_mod_full <- lme4::lmer(b ~ land_5 * region * anom_5 + (1| year), data =drop_na(df, land_5, anom_5))lag_mod <- lme4::lmer(b ~ land_5 * region + anom_5 + (1| year), data =drop_na(df, land_5, anom_5))lag_mod_sst <- lme4::lmer(b ~ anom_5 * region + (1| year), data =drop_na(df, land_5, anom_5))lag_mod_landings <- lme4::lmer(b ~ land_5 * region + (1| year), data =drop_na(df, land_5, anom_5))# Summarytidy(lag_mod) %>%gt()```### Significance```{r}# drop1(lag_mod_full, test = "Chisq")drop1(lag_mod, test ="Chisq")# drop1(lag_mod_sst, test = "Chisq")# drop1(lag_mod_landings, test = "Chisq")```### Fit & Predictions```{r}qqnorm(resid(lag_mod))qqline(resid(lag_mod))ggplot(df, aes(x = yr_num, y = b, color = region)) +geom_point(alpha =0.5) +theme_classic() +scale_color_gmri() +geom_line(data =cbind(drop_na(df, land_5, anom_5), pred =predict(lag_mod)), aes(y = pred), size =1) +# adding predicted line from mixed model theme(legend.position ="none",panel.spacing =unit(2, "lines")) # adding space between panels```### Collinearity```{r}# multicollinearity check# col_check <- check_collinearity(lag_mod_full)col_check <-check_collinearity(lag_mod)# col_check <- check_collinearity(lag_mod_landings)# col_check <- check_collinearity(lag_mod_sst)plot(col_check)```### SST Effect```{r}# SST Anomaliezsplot(ggpredict(lag_mod, ~anom_5~region), add.data = T, jitter =0.1)```### Landings Effect```{r}plot(ggpredict(lag_mod, ~land_5+region), add.data = T, jitter =0.1)```---# Median Body Weight ModelsHow has median weight changed## Contemporary SST/Landings Model```{r}# new modelnew_mod <- lme4::lmer(med_wt ~ landings * region + sst_anom + (1| year), data =drop_na(df, landings, sst_anom))# Summarytidy(new_mod) %>%gt()```### Fit & Predictions```{r}qqnorm(resid(new_mod))qqline(resid(new_mod))ggplot(df, aes(x = yr_num, y = med_wt, color = region)) +geom_point(alpha =0.5) +theme_classic() +scale_color_gmri() +geom_line(data =cbind(drop_na(df, landings, sst_anom), pred =predict(new_mod)), aes(y = pred), size =1) +# adding predicted line from mixed model theme(legend.position ="none",panel.spacing =unit(2, "lines")) # adding space between panels```### Significance```{r}drop1(new_mod, test ="Chisq")```### Collinearity```{r}# multicollinearity checkcol_check <-check_collinearity(new_mod)plot(col_check)```### SST Effect```{r}# SST Anomaliezsplot(ggpredict(new_mod, ~sst_anom+region), add.data = T, jitter =0.1)```### Landings Effect```{r}plot(ggpredict(new_mod, ~landings+region), add.data = T, jitter =0.1)```## 5-Year Smoothed Predictor Model```{r}lag_mod <- lme4::lmer(med_wt ~ land_5 * region + anom_5 + (1| year), data =drop_na(df, land_5, anom_5))# Summarytidy(lag_mod) %>%gt()```### Fit & Predictions```{r}qqnorm(resid(lag_mod))qqline(resid(lag_mod))ggplot(df, aes(x = yr_num, y = med_wt, color = region)) +geom_point(alpha =0.5) +theme_classic() +scale_color_gmri() +geom_line(data =cbind(drop_na(df, land_5, anom_5), pred =predict(lag_mod)), aes(y = pred), size =1) +# adding predicted line from mixed model theme(legend.position ="none",panel.spacing =unit(2, "lines")) # adding space between panels```### Significance```{r}drop1(lag_mod, test ="Chisq")```### Collinearity```{r}# multicollinearity checkcol_check <-check_collinearity(lag_mod)plot(col_check)```### SST Effect```{r}# SST Anomaliezsplot(ggpredict(lag_mod, ~anom_5~region), add.data = T, jitter =0.1)```### Landings Effect```{r}plot(ggpredict(lag_mod, ~land_5+region), add.data = T, jitter =0.1)```---# Region Random Effects# What if we swap Random effectsI don't really care about regions, so what if we used region as a random effect. It should also eat variability similar to year as a random effect and maybe help resolve the temperature or landings effects.Just to see:## Contemporary SST/Landings Model```{r}# new modelnew_mod <- lme4::lmer(med_wt ~ landings * sst_anom + (1| survey_area), data =drop_na(df, landings, sst_anom))# # Autocorrelative Error - Performs badly here# Ben bolker to thatnk, of course:# https://bbolker.github.io/mixedmodels-misc/notes/corr_braindump.htmlnew_mod <- nlme::lme(med_wt ~ landings * landings + sst_anom,random =~1| survey_area,data =drop_na(df, landings, sst_anom),correlation =corAR1(), method ="ML")# Summarytidy(new_mod) %>%gt()```### Fit & Predictions```{r}qqnorm(resid(new_mod))qqline(resid(new_mod))ggplot(df, aes(x = yr_num, y = med_wt, color = region)) +geom_point(alpha =0.5) +theme_classic() +scale_color_gmri() +geom_line(data =cbind(drop_na(df, landings, sst_anom), pred =predict(new_mod)), aes(y = pred), size =1) +# adding predicted line from mixed model theme(legend.position ="none",panel.spacing =unit(2, "lines")) # adding space between panels```### Significance```{r}drop1(new_mod, test ="Chisq")```### Collinearity```{r}# multicollinearity checkcol_check <-check_collinearity(new_mod)plot(col_check)```### SST Effect```{r}# SST Anomaliezsplot(ggpredict(new_mod, ~sst_anom), add.data = T, jitter =0.1)```### Landings Effect```{r}plot(ggpredict(new_mod, ~landings), add.data = T, jitter =0.1)```## 5-Year Smoothed Predictor Model```{r}# Both sst & landings lag_mod <- lme4::lmer(b ~ land_5 * anom_5 + (1| survey_area), data =drop_na(df, land_5, anom_5))# Summarytidy(lag_mod) %>%gt()```### Significance```{r}drop1(lag_mod, test ="Chisq")```### Fit & Predictions```{r}qqnorm(resid(lag_mod))qqline(resid(lag_mod))ggplot(df, aes(x = yr_num, y = b, color = region)) +geom_point(alpha =0.5) +theme_classic() +scale_color_gmri() +geom_line(data =cbind(drop_na(df, land_5, anom_5), pred =predict(lag_mod)), aes(y = pred), size =1) +# adding predicted line from mixed model theme(legend.position ="none",panel.spacing =unit(2, "lines")) # adding space between panels```### Collinearity```{r}# multicollinearity checkcol_check <-check_collinearity(lag_mod)plot(col_check)```### SST Effect```{r}# SST Anomaliezsplot(ggpredict(lag_mod, ~anom_5), add.data = T, jitter =0.1)```### Landings Effect```{r}plot(ggpredict(lag_mod, ~land_5), add.data = T, jitter =0.1)```# Variance Partitioning```{r}library(vegan)```